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The discrete-dipole approximation (DDA) is a powerful method for cal- 
culating absorption and scattering by targets that have sizes smaller 
than or comparable to the wavelength of the incident radiation. The 
, DDA can be extended to targets that are singly- or doubly-periodic. We 

I generalize the scattering amplitude matrix and the 4x4 Mueller matrix 

Qh' to describe scattering by singly- and doubly-periodic targets, and show 

. how these matrices can be calculated using the DDA. The accuracy of 

DDA calculations using the open-source code DDSCAT is demonstrated 
by comparison to exact results for infinite cylinders and infinite slabs. 
A method for using the DDA solution to obtain fields within and near 
Q^l the target is presented, with results shown for infinite slabs. 
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^ I 1. Introduction 

00 ' 

^ I Electromagnetic scattering is used to study isolated particles, but increasingly to characterize 

■ extended targets ranging from nanostructure arrays in laboratories, to planetary and asteroidal 

, regoliths. To model the absorption and scattering, Maxwell's equations must be solved for the 

^ I target geometry. 

' For scattering by isolated particles with complex geometry, a number of different theoretical ap- 

O ■ proaches have been used, including the discrete dipole approximation (DDA) [l|, Q, 0] , also known 

^ I as the coupled dipole approximation or coupled dipole method. The DDA can treat inhomogeneous 

targets and anisotropic materials, and has been extended to treat targets near substrates [1, 0]. 
Other techniques have also been employed, including the finite difference time domain (FDTD) 
5^ ; method 0,1]. 

For illumination by monochromatic plane waves, the DDA can be extended to targets that are 
spatially periodic and (formally) infinite in extent. This could apply, for example, to a periodic 
array of nanostructures in a laboratory setting, or it might be used to approximate a regolith by 
a periodic array of "target unit cells" with complex structure within each unit cell. 

Generalization of the DDA (= coupled dipole method) to periodic structures was first presented 
by Markel 0] for a 1-dimensional chain of dipoles, and more generally by Chaumet et al. 10], who 



calculated the electric field near a 2-dimensional array of parallelepipeds illuminated by a plane 



wave. Chaumet and Sentenac Hi further extended the DDA to treat periodic structures with a 
finite number of defects. 

From a computational standpoint, solving Maxwell's equations for periodic structures is only 
slightly more difficult than calculating the scattering properties of a single "unit cell" from the 
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structure. Since it is now feasible to treat targets with N > 10^ dipoles (target volume > 200 A^, 
where A is the wavelength), it becomes possible to treat extended objects with complex substruc- 
ture. 

The objective of the present paper is to present the theory of the DDA applied to scattering 
and absorption by structures that are periodic in one or two spatial dimensions. We also generalize 
the standard formalism for describing the far-field scattering properties of finite targets (the 2x2 
scattering amplitude matrix, and 4x4 Mueller matrix) to describe scattering by periodic targets. 
We show how to calculate the Mueller matrix to describe scattering of arbitrarily-polarized radi- 
ation. The theoretical treatment developed here has been implemented in the open-source code 
DDSCAT 7 (see Appendix A). 

The theory of the DDA for periodic targets is reviewed in ^ and the formalism for describing the 
far-field scattering properties of periodic targets is presented in f|3][5j Transmission and reflection 
coefficients for targets that are periodic in two dimensions are obtained in f|6l 

The applicability and accuracy of the DDA method are discussed in ^ ^7|8l where we show 
scattering properties calculated using DDSCAT 7 for two geometries for which exact solutions are 
available for comparison: (1) an infinite cylinder and (2) an infinite slab of finite thickness. The 
numerical comparisons demonstrate that, for given A, the DDA converges to the exact solution as 
the interdipole spacing d ^ 0. 

2. DDA for Periodic Targets 

The discrete-dipole approximation (DDA) is a general technique for calculating scattering and 
absorption of electromagnetic radiation by targets with arbitrary geometry. The basic theory of the 
DDA has been presented elsewhere Conceptually, the DDA consists of approximating the target 
of interest by an array of polarizable points, with specified polarizabilities. Once the polarizabilities 
are specified, Maxwell's equations can be solved accurately for the dipole array. When applied to 
finite targets, the DDA is limited by the number of dipoles N for which computations are feasible 
- the limitations may arise from large memory requirements, or the large amount of computing 
that may be required to find a solution when N is large. In practice, modern desktop computers 
are capable of solving the DDA equations, as implemented in DDSCAT 1,112, for N as large as 
~ 10^. 

Developed originally to study scattering from isolated, finite structures such as dust grains [ij, 
the DDA can be extended to treat singly- or doubly-periodic structures. Consider a collection of 

polarizable points, defining a "target unit cell" (TUC). Now consider a target consisting of a 
1-dimensional or 2-dimensional periodic array of identical TUCs, as illustrated in Figs. [T] and El 
we will refer to these as 1-d or 2-d targets, although the constituent TUC may have an arbitrary 
3-dimensional shape. For a monochromatic incident plane wave 

Einc(r, t) = Eo exp {iko ■ r - imt) , (1) 

the polarizations of the dipoles in the target will oscillate coherently. Each dipole will be affected 
by the incident wave plus the electric field generated by all of the other point dipoles. 

Let index j = 1,...,N run over the dipoles in a single TUC, and let indices m, n run over 
replicas of the TUC. The (m, n) replica of dipole j is located at 

rjmn = TjOO + "iLn + nL^ , (2) 

where Jju and L^, are the lattice vectors for the array. For 1-d targets we let m vary, but set n = 0. 
For 2-d targets, the area per TUC is 

^TUC = |Lu X L^l = LuLj; sin 9uy (3) 
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Fig. 1. (a) Target consisting of a 1-d array of TUCs and (b) showing how an infinite cylinder can be 
constructed from disk-like TUCs (lower). The M = scattering cone, with = ao, is illustrated. 



where 9uv is the angle between and L^. 

The replica dipole polarization Pjmn{t) is phase-shifted relative to Pjoo(i)- 

Pjmn(i) = Pjoo(i) exp [i(mko • L„ + nko • L^)] (4) 

Define a matrix A such that —Aj^kmn^kmn gives the electric field E at rjoo produced by an 
oscillating dipole Pkmn located at r^mn- Expressions for the 3x3 tensor elements of A have been 
presented elsewhere [e.g., A depends on the target geometry and wavelength of the incident 
radiation, but not on the target composition or on the direction or polarization state of the incident 
wave. 

Using eq. ^ we may construct a matrix A such that, for j ^ k, — A^ ^Pfcoo gives the electric 
field at Tjoo produced by a dipole Pfcoo and all of its replica dipoles Pkmm and for j = A; it gives 
the electric field at rjoo produced only by the replica dipoles: 

^j,k = (1 - <5jfc'5mO'^no) Aj^fcmnexp [i(mko • L„ nko • L^)] . (5) 



m=— oo n=— Tin 



where nmax = for 1-d targets and rimax = oo for 2-d targets, and 5ij is the Kronecker delta. For 
|m|, |n| oo, location jOO is in the radiation zone of dipole kmn, and the electric field falls off in 
magnitude only as 1/r. The sums in ([5]) would be divergent were it not for the oscillating phases of 
the terms, which ensure convergence. Evaluation of these sums can be computationally-demanding 



when koLy or fegLz are small. Chaumet et al. [lOl ] have discussed methods for efficient evaluation 
of these sums. 

We evaluate ^ numerically by introducing a factor exp[— (7A;or)^] to smoothly suppress the 
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Fig. 2. (a) Target consisting of a 2-d array of TUCs and (b) showing how an infinite slab is created 
from TUCs consisting of a single "line" of dipoles. 



contributions from large r, and truncating the sums: 



m.n 



i(mko • Lu + nko • L^) - {'ykorj^kmn)^ 



(6) 



where rj^kmn = l^^kmn — fjool and the summation is over {m,n) with rj^kmn < 2/7^0, i-e., out to 
distances where the suppression factor exp [— ~ e~^^. For given ko, L„, L^, the Aj^k 
depend only on rjoo — T^koo, and therefore only 0{8N) distinct Aj^k require evaluation. 

Ideally, one would use a very small value for the interaction cutoff parameter 7, but the number 
of terms [oc 7"^ for 1-d, or oc 7"^ for 2-d] in eq. ^ diverges as 7 — > 0. We show that 7 ~ 0.001 
ensures accurate results for the cases studied here. 

The polarizations Pjoo of the dipoles in the TUC must satisfy the system of equations 



jOO 



Einc(rj) 



(7) 



If there are dipoles in one TUC, then ([7]) is a system of 3A^ linear equations where the polariz^ 



qua i 

B 



IJ]. After A has been 



ability tensors aj are obtained from lattice dispersion relation theory 
calculated, equations ([7|) can be solved for Pjoo using iterative techniques when N ^ 1. 
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3. In the Radiation Zone 



In the radiation zone kr ^ 1, the electric field due to dipole jmn is 



E 



/cq exp (i/co|r — r 



jmn I 



'jmn 



r — r 



'■jmn 



)(r 



'-jmn) 



r — r 



jmn I 



jmn I 



r 2r • Tjmn + ^jmn 



'-jmn I 



1/2 



(8) 
(9) 



r <^ 1 



r ■ r 



jmn 



+ 



2r2 



jmn 



r • r 



jmn 



+ 



Define the unit vector = ks/fco- We seek to sum the contribution of all the dipoles to the 
electric field propagating in direction k^. At location r = rkg, the dominant contribution will be 
from dipoles located within the Fresnel zone [see, e.g., ref. Il5|, which will have a transverse radius 
Rp K, (r/Zco)^^^ • For dipoles within the Fresnel zone. 



^1 

^'jmn I 

Thus, in the radiation zone. 



(r 



'■jmn 



)(r-r 



jmn ) 



'■jmn I 



1 r 



jmn 



1 kcke 



jmn 



^3 ^ ^ 

E(r) = — ^ exp {ikQv) 1 - k^k^ ^ Pjoo X! (^^jmn 

'^Q'^ '- m.n 



jmn 



k 

mko ■ L„ + nko • L^, - k^ • rjmn + — 

2r 

-ks ■ r,oo + "T'(ko - k^) • L„ + n(ko - k^) • L„ 



jmn 



+ 



1 



2kor 



m (/cq - kg^)L^ + n [k^ - kg^)L^ + 2mn{kQLu-'Lv-ksuksvLuLy] 



+0 



(10) 

(11) 

(12) 
,(13) 



where ksu = k^ • L„/Lu, kgv = kg • hy/Ly, and terms of order (mL/r) may be neglected because 
niL/r ~ Rp/f oc r"^/^ as r — > oo. Thus, for r — > oo, the electric field produced by the oscillating 
dipoles is 



E, 



A-2 



exp (ikor) 1 - k^k^J ^ Pjoo exp (-ik^ • rjoo) ^(r, k^ 
G(r, ks) = ^ exp (i^mn) 



m,n 



m(ko - ks) • Jju + "-(ko - k 



+- 



1 



'^^(^o ~ ^su)^u + ^'^(^0 ~ ^sv)-^v ~l" 2mn[kQjju ■ Lti ~ kguksvLuLy , 



s I 



,2^7,2 



2A;or 

It is convenient to define 



AT 



FTUc(ks 



1 - ksksj X! Pjoo exp (iwt - ik^ ■ rjoo) , 

i=i 

so that the electric field produced by the dipoles is 

exp {iks ■ r — iujt) . 



Es 



kor 



-FTuc(ks)G'(r,ks 



(14) 
(15) 

(16) 
(17) 

(18) 



Ftuc depends upon the scattering direction kg, and also upon the direction of incidence ko and 
polarization Eq of the incident wave. 
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A. Isolated Finite Target: u = 

We will refer to finite isolated targets - consisting of only the dipoles in a single TUC - as targets 
that are periodic in = dimensions. For this case, we simply set G = 1 in eq. (114^ : the scattered 
electric field in the radiation zone is 



E, 



exp (ik^r — iujt) 



kor 



FTUc(ks) 



The time-averaged scattered intensity is 

. CIE.J2 



■iTT 



' TUC I 



and the differential scattering cross section is 



1 IF 



TUC I 



kl lEnP 



(19) 



(20) 



(21) 



B. Target Periodic in One Dimension: v = \ 



Without loss of generality, we may assume that targets with 1-d periodicity repeat in the y direc- 
tion. It is easy to see from eq. (|15p and (jl6p that G = except for scattering directions satisfying 



k 



sy 



koy + M— 



M = 0,±1,±2,..., 



where energy conservation 



k1 = k^) limits the allowed values of the integer M: 
i-koy - ko)^ < M < i-koy + ko)^ . 



(22) 



(23) 



If (/cq + \koy\)Ly < 2iT, then only M = scattering is allowed. Define polar angles ao and Og for 
the incident and scattered radiation, so that 



koy 



^sy 



ko cos ao 
ko cos as 



(24) 
(25) 



For each allowed value of kgy, the scattering directions define a cone: 



ksyf + {ki - kiy/^ 



sma. 



sm ao 



(ko — y cos ao) cos C + y x ko sin ( 



(26) 



where C is an azimuthal angle measured around the target axis y- The sum for G(r, k^) is [since 
M must be an integer - see eq. (I22p ] 



G= ^ exp(i<I> 



mOj 



oo 

J2 



m=—oo 



m=—oo 

lim 



-2TTiMm + i 



m 

2kor 



(/Cg kgy)Ly 



dm exp 



i{l + ie) 2(1,2 1.2 \t2 



{2mkor)^/^ {2'Kikor)'^/^ 



(^0 ~ ^lyY^'^^y ^oLy sin 



and the scattered electric field 



/ 27r A exp {ik^r - iivt) ^ 
E« = ttt: 1. T — FTUc(k 



\korJ 



koLy sin a^ 



(27) 
(28) 
(29) 

(30) 
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shows the expected r~^/^ behavior far from the scatterer (the distance from the cyhnder axis is 
R = r sina^). 

For each allowed value of k^, the total time-averaged scattered power Psca., per unit length of 
the target, per unit azimuthal angle (, may be written 



dLdC Sir dLdC 
where the differential scattering cross section is 

d^Csca Svr Id^Psc^ Svr lEpc 



R sin Os (32) 



dLdC |Eop c dLdC |Eopc Svr 

(33) 



2tt [FtucI^ 



kiLl |Eo|2 ■ 
C. Target Periodic in Two Dimensions: v = 2 

For targets that are periodic in two dimensions, it is apparent from eq. (jl5p and (jl6p that G = 
unless 

(k, -ko)-L„ = 2tiM , M = 0,±l,±2,.., (34) 
(k5-ko)-L^ = 2ttN , iV = 0,±1,±2, .... (35) 

The 2-d target constitutes a diffraction grating, with scattering allowed only in directions given by 
(j34l35p . It is convenient to define the reciprocal lattice vectors 

27rx X 27rx x L„ 

u = ^— — — , V = ^— — — (36) 

X • (L„ X L^) X • (L^, X L„) 

The wave vector transverse to the surface normal is 

k.± = ko± + Mu + iVv . (37) 

Energy conservation requires that 

kl^ = kl- |ko_L + Mu + iVvp > . (38) 

For any {M,N) allowed by eq. ([38]) . there are two allowed values of ksx, differing by a sign; one 
(with ksxkox > 0) corresponds to the {M,N) component of the transmitted wave, and the other 
(with ksxkox < 0) to the {M,N) component of the reflected wave. Define 

sinag = -^^^ , (39) 
ko 

sin a, = M . (40) 

Kq 

Note that ao = "^f^ for normal incidence, and oq ^ for grazing incidence. For k<j satisfying eq. 
(I34I35I) . we have 

G = ^exp(i$m„) 

m,n 

= ^ exp I {kl - kl^ Liw? + {kl - kl^ Lln^ + 2 (^k^Lu ■ 'Lv - ksuksvLuL^^ mn 
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lim 



X exp 



lim 



dm 

oo 

i(l + ie 



dn 



2kor 



1 



e-^0+ AtUC 

2Tiir 



'"suJ-^u 
oo 



dy I dzexp 



oo 



i{l+ie) 



2fcor 



koA-cuc sin 
The scattered electric field is 



E, 



27ri exp (ik^ 



fco^TUC smus 

Note that |Es| is independent of distance x = r sina^ from the target, as expected for a target that 
is infinite in two directions. For pure forward scattering, = ko, we must sum the incident wave 
Einr and the radiated wave E,: 



' TUC 



(41) 



(42) 



E = exp (iko • r — iut) 



Eo + 



27rzFTuc(ks = ko) 



^o^TUC sin Us 

The cross section per unit target area A for scattering into direction {M,N) is (for k^ ^ ko): 
dC7sca(M,iV) lEpsina, 



dA 



|Eop siuQo 
47r2 



(43) 



(44) 



'TUC I 



(45) 



where Cgca can be evaluated for either transmitted or reflected waves. For the special case k^ = kg, 
the transmission coefficient T(M,N) is obtained from the total forward-propagating wave (|43p : 

2 



T(0,0) 



1 



|EoP 



Eo + 



27ri FTUc(ks=k 



^o^TUC smao 

4. Scattering Amplitude Matrices 5^-'^'^^ 
A. Isolated Finite Targets: v = Q 

In the radiation zone, the scattered electric field is related to the incident electric field via the 2x2 



(46) 



scattering amplitude matrix [16'], defined so that 

I ■ e^ii \ 
I E^ • e^_L / 



iexp(ik, -r-io;*) ( sf^^ sf'^^ 

(Od) „{0c() 



where 



kor 



ko X k. 







\ks X kol 
ko X eij_ 



1 



ks X ko 
- (k. • ko)2 



Eo -ei 
Eo • ei_ 



(47) 

(48) 
(49) 

(50) 

l-(k.-ko)2 

are the usual conventions for the incident and scattered polarization vectors parallel and perpen- 
dicular to the scattering plane [see, e.g., §3.2 of ref. Il6l|. 



-II 



= kc X e 



ks - (kg • ko)ko 

1 - (ks • ko)2 
-ko + {ks ■ ko)ks 
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B. Target Periodic in One Dimension: v = 1 

For targets with 1-d periodicity, it is natural to generalize the scattering amplitude matrix, so 
that - for directions for which scattering is allowed - the scattered electric field at a distance 
R = r sin as from the target is 



Es • egii \ zexp(iks • r — iwt) / So^^'^ S'o^'^^ A / Eq • 



(51) 



for ks satisfying eq. (|22ll26p . 

C. Target Periodic in Two Dimensions: v = 2 

For targets with 2-d periodicity, it is natural to generalize the scattering amplitude matrix so that, 
for directions k^ 7^ ko for which scattering is allowed, we write 



for ks satisfying eq. (|37ll38p . 

For the special case of forward scattering (M = = and k^ = kg), where the scattering plane 
is not simply defined by kg and k^, it is natural to use kg and the target normal x to define the 
scattering plane. Thus 



ko X ksj 



|ko X ks_L| 

with ejii and e^y defined by (|49l50p . For r ^ 00 



(53) 



D. Far-Field Scattering Amplitude Matrices 

The scattei 
V = 0,1,2: 



The scattering amplitude matrices s\^'^'' (k-s) are directly related to the Ftuc for the three cases, 



S^{^^ = C,esi_ ■ FTuc(k., Eo = e,^) (55) 

Si"''^ = ae,|| • FTuc(k„Eo = e,,,) (56) 

Si"''^ = ae,|| • FTuc(k., Eo = e,^) (57) 

= C,esi_ ■ FTuc(k., Eo = e,,,) (58) 



Co = -i (59) 
C,^-(^)'"-l- (60) 

= WT^- (61) 



for finite targets (Co), and targets that are periodic in one or two dimensions (Ci or C2) 
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5. Far-Field Scattering Matrix for Stokes Vectors 

For a given scattering direction k<j, the 2x2 complex amplitude matrix S^'^'^\]<.s) fully characterizes 
the far-field scattering properties of the target. The far-field scattering properties of an isolated 
finite target are characterized by the 4x4 dimensionless Mueller matrix S^^\ with the Stokes 
vector of radiation scattered into direction k^ at a distance r from the target given by 

1 (Od) 
/3 1 

where /inc,/3 = {Ii Qi U, F)inc is the Stokes vector for the radiation incident on the target. For 1-d 
targets, we define the dimensionless scattering matrix S^^^^ by 

1 (Id) 

Isca.,a = 2^ hncfi , (63) 

where R is the distance from the one-dimensional target. For 2-d targets, we define 5*^^^ by 

4 

-^sca.a = X! hncfi ■ (64) 

(3=1 

The 4x4 scattering intensity matrix S^^^^ is obtained from the scattering amplitude matrix elements 

gb^d) _ Except for the special case of forward scattering (k^ = kg) for 2-d targets, the equations are 
the same as eq. (3.16) of Bohren and Huffman [3]. For example 



s'^-> = 1 {^\sr'v + \sr'\ + i^rf + 1^4 tj (65) 

4f = \ I' - iM'"^ I' - l^r^ I' + I'^r^ I') (66) 

5if)=Im(4'"'^5("''^*-5i"'')5f''^*) . (67) 

For the special case of forward scattering (k^ = ko) for 2-d targets, it is necessary to replace S^'^^ 
and Sj^*^^ with (S'|^'^^ — i) and — i) [cf. Eq. Thus, for example, 

.{2d)/, N _ 1 (\c{'^d) .|2 , |c.{2rf) ■\2 I |c(2d)|2 , |c(2d)|2 



5}r(k«=ko) = ^ [\sr>-i\' + i5r'-ii^ + i^rr + 1'54 tj • (68) 

6. Transmission and Reflection Coefficients for 2-D Targets 

For targets with 2-d periodicity, it is natural to define generalized transmission and reflection 
coefficients for the Stokes vectors: for scattering order {M,N), /sca,a = J2i3'^af3{M, N)I[ac,f3 is the 
Stokes vector component a for radiation with ksxki^cx > 0, and Rafj[M,N) is the fraction of the 
incident Stokes vector component fi that emerges in Stokes vector component a with ksxki-ac,x < 0. 
These can be related to the S^^^^ : 

R^p{M,N) = ^-^S^^^^ foTksxkox<0 , (69) 
smao ^ 

r„^(M,iV) = ^^^^d) for A:,,A:o:. > , (70) 
sm ao ^ 
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Fig. 3. Scattering by an infinite cylinder with diameter D and m = 1.33 + O.Oli, for radiation with 
X = ttD/X = 50 and incidence angle ao = 60°. (a) S^'^\ Solid curve: exact solution. Broken 
curves: DDA results for D/d= 256, 360, and 512 {N = 51676, 102036, 206300 dipoles per TUC); 
(b) fractional error in S'j}'^^ (DDA). (c) s!^i'^'^ ; (d) error in s!^i^\ 
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Fig. 4. Scattering by an infinite cylinder with diameter D and m = 1.33 + O.Oli, for radiation with 
ttD/X = 50, and incidence angle ao = 60°. (a) Exact solution (solid curve) and DDA results for 
D/d = 512 and various values of the interaction cutoff parameter (b) fractional error in S^-^ ^5 
(c,d) same as (a,b), but expanding the region < C < 20°. For this case, results computed with 
7 = 0.002 and 0.001 are nearly indistinguishable. 
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Fig. 5. Light scattered by an infinite cylinder with m = 2 + i, for radiation with x = 2ttR/X = 25 
and incidence angle ao = 60°. (a) S^'^\ Solid curve: exact solution. Broken curves: DDA results 
for D/d= 128, 180, and 256 {N = 12972, 25600, 51676 dipoles per TUC). (b) Fractional error in 
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The fraction of the incident power that is absorbed by the target is 

P bJArea ^ ^ _ ^ ^ ^) + y ^y)] , (71) 

|Eo|^csmao/87r j^^i ^mc,i 

where line,/? is the Stokes vector of the incident radiation. 

For unpolarized incident radiation, Rn (M, A'^) is the fraction of the incident power that is 
reflected in diffraction component (M, N), Tii(M, N) is the fraction that is transmitted in compo- 
nent (M, A^), and 1 — Sm7v[-^ii(-^) ^) + ^ii(-^) ^)] is tlie fraction of the incident power that is 
absorbed. 

7. Example: Infinite Cylinder 

DDSCAT 7 has been used to calculate scattering and absorption by an infinite cylinder consisting 
of a periodic array of disks of thickness d and period Ly = d (where d is the interdipole spacing). 
Fig. [3^ shows S^f^ for refractive index m = 1.33 + O.Oli and ttD/X = 50 {D is the cylinder 
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Table 1. CPU time to calculate scattering by m = 1.33 + O.Oli infinite cylinders on single-core 2.4 
GHz AMD Opteron model 250 



ttD/X 


N 


7 


calc. A 


solution 


scat. 


Total 








(min) 


(min) 


(min) 


(min) 


25 


51676 


0.005 


3.29 


17.6 


0.65 


22.2 


25 


51676 


0.001 


16.2 


17.8 


0.65 


35.3 


50 


102036 


0.005 


4.58 


59.7 


1.27 


66.8 


50 


102036 


0.001 


22.8 


45.2 


1.09 


70.2 


50 


206300 


0.005 


13.2 


254. 


2.88 


273. 


50 


206300 


0.001 


66.0. 


292. 


2.76 


364. 



diameter and A the wavelength of the incident radiation), and incidence angle ao = 60°. Because 
/co(l + I cosao|)d < 2tt, equations (p2t \23\i allow only M = scattering, with Og = oq. Also shown 
is the exact solution, calculated using a code written by D. Mackowski (private communication). 
Light scattering by cylinders is generally described by scattering amplitudes Tj; in Appendix [Bl we 
provide expressions relating these Tj to the Si used here. Fig. \3]p shows the fractional error in S'j^'^^ 
calculated using DDSCAT. As d is decreased, the errors decrease. Excellent accuracy is obtained 
when the validity criterion [1] \m\kd ^ 0.5 is satisfied: the fractional error in is typically less 
than a few %, except near deep minima in 5ii. 

Fig. [3h shows 82^^ , characterizing scattering of unpolarized light into the Stokes parameter Q 
(S'21 < corresponds to linear polarization perpendicular to the scattering plane). DDSCAT 7 and 
the exact solution are in very good agreement when \m\kd ^ 0.5. Note that although the error in 
S2Y\6 = 0) ?a 2 is large compared to 5*21(0) = —6, this is small compared to s[Y\o) ^ 1500: the 
scattered radiation is only slightly polarized. 

The results in Fig. [3] were obtained using 7 = 0.001 to truncate the integrations. To see how 
the results depend on 7, Figure H] shows s[Y^ computed for the problem of Fig. [3] but using 
different values of 7. For azimuthal angles C > 20°, the results for 7 = 0.005 and 0.001 are nearly 
indistinguishable; the difference between the computed result and the exact solution is evidently 
due to the finite number of dipoles used, rather than the choice of cutoff parameter 7. However, 
the results for forward scattering are more sensitive to the choice of 7, as is seen in Fig. [31:, d: it is 
necessary to reduce 7 to 0.001 to attain high accuracy in the forward scattering directions. 

Table [1] gives the CPU times to calculate A, to then iteratively solve the scattering problem to 
a fractional error < 10~^ (using double-precision arithmetic), and finally to evaluate the scattering 
intensities, for several of the cases shown in Figs. [3] and [H For most cases the CPU time is 
dominated by the iterative solution using the conjugate gradient algorithm. While the time required 
to evaluate A might be reduced using the strategies suggested by [13], this step is generally a 
subdominant part of the computation for targets with kd^ 0.1. 

The above results have been for a weakly-absorbing cylinder. To confirm that the DDA can 
be applied to strongly-absorbing material. Fig. [5] shows scattering calculated for a cylinder with 
m = 2 + i and x = ttD/X = 25. Once again, the accuracy is very good, with small fractional errors 
provided \m\kd ^ 0.5. 

8. Example: Plane-Parallel Slab 

Consider a homogeneous plane-parallel slab with thickness h and refractive index m. Radiation 
incident on it at angle of incidence 9i will either be specularly refiected or transmitted. The 
refiection and transmission coefficients R and T can be calculated analytically, taking into account 
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Fig. 6. Transmission and reflection coefficients for radiation with wavelength A incident at angle 
9i = (7r/2 — ao) = 40° relative to the normal on a slab with thickness h, incident E || and _L to the 
scattering plane, as a function of Re(m)/i/A. (a) Nonabsorbing slab with m = 1.5 (b) Absorbing 
slab with m = 1.5 + 0.02i. Solid curve: exact solution. Symbols: results calculated with the DDA 
using dipole spacing d = h/W, h/20, and h/AO. 



multiple reflections within the slab jl5(]. With an exact solution in hand, we can evaluate the 
accuracy of the DDA applied to this problem. Figure [6] shows results for two cases: a dielectric 
slab with m = 1.50, and an absorbing slab, with m = 1.50 + 0.021 

DDSCAT 7 was used to calculate reflection, transmission, and absorption by an inflnite slab, 
generated from a TUC consisting of a single line of dipoles extending in the x direction, with Ly = 
Lz = d. The selection rules ()34|35|38p allow only M = = 0: transmission or specular reflection. 
The reflection and transmission coefficients for radiation polarized parallel or perpendicular to the 
plane containing k and the surface normal are 

-ko.) + s[f{ks^ = ko,) (72) 
-ko.) - Sg^'^ks. = ko,) (73) 
ko.) + S[l''\ks^ = ko^) (74) 

kox) - S[l'''\ksx = kox) ■ (75) 

The DDA results are in excellent agreement with the exact results when the validity condition 
\m\kQd < 0.5 is satisfied, but results with moderate accuracy are obtained even when \m\kQd ~ 1. 

9. Near-Field Evaluation 

The polarizations Pjoo can be used to calculate the electric and magnetic fields at any point, 
including within or near the target, using the exact expression for E and B from a point dipole, 
modified by a function (p: 

E(r, t) = e ^ ' ^ ' rjz ^ (pi-^jmrOl k^TLj^nn X {Fjmn ^ R-jmn) 

m,n \^jmn\ L 



-^11 — "^11 {ksx — 

= ^ [ksx — 

rp _ qi^d.),, _ 

II — 11 \"'sx — 

T± = {kgx = 
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Fig. 7. |E^|/|Eq| along two tracks normal to the dielectric slab of Fig. [6)3, for slab thickness 
h = 0.2A, incidence angle Oj = 40°, and incident polarizations || and _L to the scattering plane (see 
text). Results were calculated using eq. ([76]) with the slab represented by = 10 and A^^^ = 20 
dipole layers (i.e., dipole spacing d = O.l/i and 0.05/i). The circles along track 1 are at points where 
dipoles are located. 



+ - 



[1 — iknR. 



B(r,t) 



-iuit 



tj2 



jmn J 



Rjmn^jmn 



+ Eq exp(zko • r — iujt) (76) 



2 eyip{ikoRjr^ 



-(f>{Rjmn) (R-jmn ^ Pj'mji) I 1 



J m,n 

+ ko X Eq exp(iko • r — iujt) 

f 1 for R>d 

^ 1 (R/d)^ for R<d 



iknR 



exp 



(77) 
(78) 

(79) 



The function (f){R) smoothly suppresses the (oscillating) contribution from distant dipoles in order 
to allow the summations to be truncated, just as in eq. ^ for evaluation of Aj fc. If r is within 
the target or near the target surface, the summations over (m, n) are limited to \Rjmn\ < ^/"jkQ. 
The (R/d)^ factor suppresses the R~^ divergence of E as r approaches the locations of individual 
dipoles, and at the dipole locations results in E that is exactly equal to the field that is polarizing the 
dipoles in the DDA formulation. Evaluation of eq. (j76|. [77|) is computationally-intensive, because 
the summations EjE^„ typically have many terms. 

To illustrate the accuracy, we consider the infinite slab of Fig. [6h), with refractive index m = 
1.5-|-0.02i and radiation incident at an angle 6i = 40°. Figure [7] shows the time-averaged |Ep/|EoP 
for slab thickness h = 0.2A - near a minimum in transmission, and a maximum in reflection (see 
Fig. [6)3. The program DDfield (see Appendix [A]) was used to evaluate E along two lines normal to 
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the slab: track 1 passes directly through dipole sites, and track 2 passes midway between the four 
nearest dipoles as it crosses each dipole layer. The E fields calculated along tracks 1 and 2 are very 
similar, although of course not identical. Within the slab, |E| along track 2 tends to be slightly 
smaller than along track 1, but for this example the difference is typically less than Figure 
[7]shows results for the slab represented by h/d = = 10 and 20 dipole layers (with \m\kd = 0.19 
and 0.094, respectively). 

Even for Nx = 10, the electric field at points more than a distance d from the edge is obtained 
to within ~ 2% accuracy at worst, which is perhaps not surprising because, as seen in Figure 
El the calculated transmission and reflection coefficients are very accurate. The discontinuity in 
\E\'^ at the boundary is spread out over a distance ~ d. The DDA obviously cannot reproduce 
field structure near the target surface on scales smaller than the dipole separation d, but fields on 
scales larger than d appear to be quite accurate. DDSCAT and DDfield should be useful tools for 



studying electromagnetic fields around arrays of nanostructures, such as gold nanodisks jlTI. Il8l|. 
10. Summary 

The principal results of this study are as follows: 

1. The DDA is generalized to treat targets that are periodic in one or two spatial dimensions. 
Scattering and absorption of monochromatic plane waves can be calculated using algorithms 
that parallel those used for finite targets. 

2. A general formalism is presented for description of far- field scattering by targets that are 
periodic in one or two dimensions using scattering amplitude matrices and Mueller matrices 
that are similar in form to those for finite targets. 

3. The accuracy of the DDA for periodic targets is tested for two examples: infinite cylinders 
and infinite slabs. The DDA, as implemented in DDSCAT 7, is accurate provided the validity 
criterion |m|A:d ^ 0.5 is satisfied. 

4. We show how the DDA solution can be used to evaluate E and B within and near the target, 
with calculations for an infinite slab used to illustrate the accuracy of near-field calculations. 
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Appendix A: DDSCAT 7 and DDfield 

The theoretical developments reported here have been implemented in a new version of the open- 
source code DDSCAT (http: //www. astro. princeton.edu/'-^draine/DDSCAT. html). DDSCAT 7 
is written in Fortran 90, with dynamic memory allocation and the option to use either single- or 
double-precision arithmetic. DDSCAT 7 includes options for various target geometries, including 
a number of periodic structures. A program DDfield for near-field calculations is also provided. 

DDSCAT 7 offers the option of using an implementation of BiCGstab with enhancement to 
maintain convergence in finite precision arithmetic [l^ . The matrix- vector multiplications AP are 
accomplished efficiently using FFTs Documentation for DDSCAT is available from ArXiv 



2ll |. with additional information available from http://ddscat.wikidot.com. 
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In addition to differential scattering cross sections, DDSCAT reports dimensionless "efficiency 
factors" Qx = Ca;(TUC)/7ra^g for scattering and absorption, where Ca;(TUC) is the total cross 
section for scattering or absorption per TUC, normalized by vrOgg, where aes = (3VTUc/4vr)^/^ is 
the radius of a sphere with volume equal to the solid volume Vtuc in one TUC. 

In the case of one-dimensional targets, with periodicity Ly in the y direction, the absorption, 
scattering, and extinction cross sections per unit target length are 

dCx 1 „ 2 



dL Ly 



QxTraiff (Al) 



for X = ext, sea, and abs, where Qx are the efficiency factors calculated by DDSCAT. 

In the case of two-dimensional targets, with periodicities L„ and Ly, the absorption, scattering, 
and extinction cross sections per unit target area are 

dCx ^ QxTTajg ^^2) 
dA LuLy sm Ouv 
Appendix B: Relation Between Si and Tj for Infinite Cylinders 

The analytic solution for infinite cylinders decomposes the incident and scattered radiation into 
components polarized parallel and perpendicular to planes containing the cylinder axis and the 
propagation vector ko or kc,. These polarization basis states differ from the choice that is usual for 
scattering by finite particles, where it is customary to decompose the incident and scattered waves 
into components polarized parallel and perpendicular to the scattering plane - the plane containing 
ko and k^. 

In the notation of Bohren and Huffman [l^ , the radiation scattered by an infinite cylinder can 
be written 



E. -e 



(ck) 




iexp(ik,.r-it^t) — — — I ^ „ :L (B1) 



where R is the distance from the cylinder axis, a is the angle between kg and the cylinder axis 
c, and superscript (ck) denotes polarization vectors parallel or perpendicular to planes containing 
the cylinder axis c and either ko or k<j. The azimuthal angle ^ is measured around the cylinder 
axis c, with C = for forward scattering. The scattering angle 6 = arccos[ko • k^] is 

9 = arccos 1 — (1 — cos Q sin^ a . 

The scattering amplitude matrix elements Tj appearing in ()B1|1 can be related to the matrix 
elements Si appearing in eq. ([5T|) : 

2i 



vr sm a 



ATB-' , (B3) 



g(ld) gild) 
(Id) cild) 



S\ Si 




e 



1 , (B5) 



i'^^) '• '-(ck) \ -I / / n\ 

" " ' ^s± 1 _ 1 ( — cot a(l — cos ^) sinasin(, 

.(ck) ^ -(cA;) 



• e^ll e^x • e^^ / sint' y — sinasin^ — cota(l — cos 

-> (ck) - [ck) \ / \ 

^ I -i|r ^iW^i± I _ 1 / cota(l-cos6') sinasinC \ ^ 

^ ^ \± ■ e^i^l'^ e,^ ■ 4? j ~ sine [ -sinasinC cota(l-cos0)y ' ^ ' 

_ sinO I cota(l — COS0) —sinasinC \ ^^^^ 

cot^ a(l — cos0)2 + sin^ asin^ C V sinasinC cota(l — cos0) / 



DDA for periodic targets 



18 



References 

1. E. M. Pur cell and C. R. Pennypacker, "Scattering and Absorption of Light by Nonsptierical 
Dielectric Grains," Astrophys. J. 186, 705-714 (1973). 

2. B. T. Draine, "The discrete-dipole approximation and its application to interstellar graphite 
grains," Astrophys. J. 333, 848-872 (1988). 

3. B. T. Draine and P. Flatau, "Discrete-dipole approximation for scattering calculations," 
J. Opt. Soc. Am. All, 1491-1499 (1994). 

4. B. T. Draine, "The Discrete Dipole Approximation for Light Scattering by Irregular Targets," 
in "Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications," , 
M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, eds. (San Diego: Academic Press, 2000), 
pp. 131-145. 

5. R. Schmchl, B. M. Nebcker, and E. D. Hirlcman, "Discrete-dipole approximation for scatter- 
ing by features on surfaces by means of a two-dimensional fast fourier transform technique," 
J. Opt. Soc. Am. A 14, 3026-3036 (1997). 

6. M. Paulus and O. J. F. Martin, "Green's tensor technique for scattering in two-dimensional 
stratified media," Phys. Rev. E 63, 066615 (2001). 

7. P. Yang and K. N. Liou, "Finite Difference Time Domain Method for Light Scattering by 
Nonspherical and Inhomogeneous Particles," in "Light Scattering by Nonspherical Particles: 
Theory, Measurements, and Applications," , M. I. Mishchenko, J. W. Hovenier, and L. D. 
Travis, eds. (San Diego: Academic Press., 2000), pp. 173 221. 

8. A. Taflove and S. G. Hagness, Advances in Computational Electrodynamics: the Finite- 
Difference Time- Domain Method (Artech House, Boston, 2005). 

9. V. A. Markel, "Coupled-dipole Approach to Scattering of Light from a One-dimensional 
Periodic Dipole Structure," Journal of Modern Optics 40, 2281-2291 (1993). 

10. P. G. Ghaumet, A. Rahmani, and G. W. Bryant, "Generalization of the coupled dipole method 
to periodic structures," Phys. Rev. B 67, 165404 (2003). 

11. P. C. Ghaumet and A. Sentenac, "Numerical simulations of the electromagnetic field scattered 
by defects in a double-periodic structure," Phys. Rev. B 72, 205437-20544 (2005). 

12. B. T. Draine and P. Flatau, "User Guide for the Discrete Dipole Approximation Gode 
DDSCAT.6.1," http://arXiv.or/abs/astro-ph/ArXiv/0409262 (2004). 

13. B. T. Draine and J. Goodman, "Beyond Glausius-Mossotti - Wave propagation on a po- 
larizable point lattice and the discrete dipole approximation," Astrophys. J. 405, 685-697 
(1993). 

14. D. Gutkowicz-Krusin and B. T. Draine, "Propagation of Electromagnetic Waves on a Rect- 
angular Lattice of Polarizable Points," http://arXiv.org/abs/astro-ph/0403082 (2004). 

15. M. Born and E. Wolf, Principles of Optics (Gambridge Univ. Press, Gambridge, 1999). 

16. G. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles 
(Wiley, New York, 1983). 

17. Z. N. Utegulov, J. M. Shaw, B. T. Draine, S. A. Kim, and W. L. Johnson, "Surfacc-plasmon 
enhancement of Brillouin light scattering from gold-nanodisk arrays on glass," in "Plasmonics: 
Metallic Nanostructures and Their Optical Properties V.", Edited by Stockman, Mark I.. 
Proceedings of the SPIE, vol. 6641 (2007), 66411M. 

18. W. L. Johnson, S. A. Kim, Z. N. Utegulov, and B. T. Draine, "Surfacc-plasmon fields in 
two-dimensional arrays of gold nanodisks," submitted for publication in SPIE 2008 Optics 
and Photonics (2008). 

19. M. A. Botchev, subroutine zbcg2, http://www.math.uu.nl/people/vorst/zbcg2.f90 

(2001). 

20. J. J. Goodman, B. T. Draine, and P. J. Flatau, "Application of fast-Fourier transform tech- 



DDA for periodic targets 



19 



niques to the discrete dipole approximation," Optics Lett. 16, 1198-1200 (1990). 
21. B. T. Draine and P. Flatau, "User Guide for the Discrete Dipole Approximation Code 
DDSCAT 7.0," http://arXiv.org/abs/astro-ph/0809.0337 (2008). 



